Quantitative Big Imaging

Kevin Mader
10 March 2016

Basic Segmentation and Discrete Binary Structures

Course Outline

  • 25th February - Introduction and Workflows
  • 3rd March - Image Enhancement (A. Kaestner)
  • 10th March - Basic Segmentation, Discrete Binary Structures
  • 17th March - Advanced Segmentation
  • 24th March - Analyzing Single Objects
  • 7th April - Analyzing Complex Objects
  • 14th April - Spatial Distribution
  • 21st April - Statistics and Reproducibility
  • 28th April - Dynamic Experiments
  • 12th May - Scaling Up / Big Data
  • 19th May - Guest Lecture - High Content Screening
  • 26th May - Guest Lecture - Machine Learning / Deep Learning and More Advanced Approaches
  • 2nd June - Project Presentations

Lesson Outline

  • Motivation
  • Qualitative Approaches
  • Thresholding
    • Other types of images
    • Selecting a good threshold
  • Implementation
  • Morphology
  • Contouring / Mask Creation

Applications

  • Simple two-phase materials (bone, cells, etc)
    • Beyond 1 channel of depth
  • Multiple phase materials
  • Filling holes in materials
  • Segmenting Fossils
  • Attempting to segment the cortex in brain imaging Cortex Image

Literature / Useful References

  • Jean Claude, Morphometry with R
  • John C. Russ, “The Image Processing Handbook”,(Boca Raton, CRC Press)
    • Available online within domain ethz.ch (or proxy.ethz.ch / public VPN)

Motivation: Why do we do imaging experiments?

  • Exploratory
    • To visually, qualitatively examine samples and differences between them
    • No prior knowledge or expectations
  • To test a hypothesis
    • Quantitative assessment coupled with statistical analysis
    • Does temperature affect bubble size?
    • Is this gene important for cell shape and thus mechanosensation in bone?
    • Does higher canal volume make bones weaker?
    • Does the granule shape affect battery life expectancy?
  • What we are looking at Standard Cell, http://en.wikipedia.org/wiki/File:Average_prokaryote_cell-_en.svg
  • What we get from the imaging modality Single Cell

To test a hypothesis

  • We perform an experiment bone to see how big the cells are inside the tissue \[ \downarrow \] Bone Measurement

2560 x 2560 x 2160 x 32 bit = 56GB / sample

  • Filtering and Preprocessing!
    \[ \downarrow \]
  • 20h of computer time later …
  • 56GB of less noisy data
  • Way too much data, we need to reduce

What did we want in the first place

Single number:

  • volume fraction,
  • cell count,
  • average cell stretch,
  • cell volume variability

Why do we perform segmentation?

  • In model-based analysis every step we peform, simple or complicated is related to an underlying model of the system we are dealing with
  • Occam's Razor is very important here : The simplest solution is usually the right one
    • Bayesian, neural networks optimized using genetic algorithms with Fuzzy logic has a much larger parameter space to explore, establish sensitivity in, and must perform much better and be tested much more thoroughly than thresholding to be justified.
    • We will cover some of these techniques in the next 2 lectures since they can be very powerful particularly with unknown data

Review: Filtering and Image Enhancement

  • This was a noise process which was added to otherwise clean imaging data
  • \[ I_{measured}(x,y) = I_{sample}(x,y) + \text{Noise}(x,y) \]
  • What would the perfect filter be
    • \[ \textit{Filter} \ast I_{sample}(x,y) = I_{sample}(x,y) \]
    • \[ \textit{Filter} \ast \text{Noise}(x,y) = 0 \]
    • \[ \textit{Filter} \ast I_{measured}(x,y) = \textit{Filter} \ast I_{real}(x,y) + \textit{Filter}\ast \text{Noise}(x,y) \rightarrow \bf I_{sample}(x,y) \]
  • What most filters end up doing \[ \textit{Filter} \ast I_{measured}(x,y) = 90\% I_{real}(x,y) + 10\% \text{Noise}(x,y) \]
  • What bad filters do \[ \textit{Filter} \ast I_{measured}(x,y) = 10\% I_{real}(x,y) + 90\% \text{Noise}(x,y) \]

Qualitative Metrics: What did people used to do?

  • What comes out of our detector / enhancement process Single Cell
  • Identify objects by eye
    • Count, describe qualitatively: “many little cilia on surface”, “long curly flaggelum”, “elongated nuclear structure”
  • Morphometrics
    • Trace the outline of the object (or sub-structures)
    • Can calculate the area by using equal-weight-paper
    • Employing the “cut-and-weigh” method

Segmentation Approaches

They match up well to the world view / perspective

Approaches

Model-Based

  • \( \rightarrow \) Experimentalist
  • Problem-driven
    • Top-down
    • Reality Model-based

Machine Learning Approach

  • \( \rightarrow \) Computer Vision / Deep Learning
  • Results-driven

Model-based Analysis

Traditional Imaging

  • Many different imaging modalities ( \( \mu \textrm{CT} \) to MRI to Confocal to Light-field to AFM).
  • Similarities in underlying equations
    • different coefficients, units, and mechanism

\[ I_{measured}(\vec{x}) = F_{system}(I_{stimulus}(\vec{x}),S_{sample}(\vec{x})) \]

Direct Imaging (simple)

  • \( F_{system}(a,b) = a*b \)
  • \( I_{stimulus} = \textrm{Beam}_{profile} \)
  • \( S_{system} = \alpha(\vec{x}) \)

\( \longrightarrow \alpha(\vec{x})=\frac{I_{measured}(\vec{x})}{\textrm{Beam}_{profile}(\vec{x})} \)

Single Cell

Nonuniform Beam-Profiles

In many setups there is un-even illumination caused by incorrectly adjusted equipment and fluctations in power and setups

Gradient Profile

Frequently there is a fall-off of the beam away from the center (as is the case of a Gaussian beam which frequently shows up for laser systems). This can make extracting detail away from the center challenging

Gaussian Beam

Absorption Imaging (X-ray, Ultrasound, Optical)

  • For absorption/attenuation imaging \( \rightarrow \) Beer-Lambert Law \[ I_{detector} = \underbrace{I_{source}}_{I_{stimulus}}\underbrace{\exp(-\alpha d)}_{S_{sample}} \]

    • Different components have a different \( \alpha \) based on the strength of the interaction between the light and the chemical / nuclear structure of the material \[ I_{sample}(x,y) = I_{source}\exp(-\alpha(x,y) d) \] \[ \alpha = f(N,Z,\sigma,\cdots) \]
  • For segmentation this model is:

    • there are 2 (or more) distinct components that make up the image
    • these components are distinguishable by their values (or vectors, colors, tensors, …)

Attenuation to Intensity

Image Histogram

Where does segmentation get us?

  • We convert a decimal value (or something even more complicated like 3 values for RGB images, a spectrum for hyperspectral imaging, or a vector / tensor in a mechanical stress field)
  • to a single, discrete value (usually true or false, but for images with phases it would be each phase, e.g. bone, air, cellular tissue)

  • 2560 x 2560 x 2160 x 32 bit = 56GB / sample \[ \downarrow \]

  • 2560 x 2560 x 2160 x 1 bit = 1.75GB / sample

Applying a threshold to an image

Start out with a simple image of a cross with added noise \[ I(x,y) = f(x,y) \]

The intensity can be described with a probability density function \[ P_f(x,y) \] Probability density function

Applying a threshold to an image

By examining the image and probability distribution function, we can deduce that the underyling model is a whitish phase that makes up the cross and the darkish background

With Threshold Overlay

Applying the threshold is a deceptively simple operation

\[ I(x,y) = \begin{cases} 1, & f(x,y)\geq0.5 \\ 0, & f(x,y)<0.5 \end{cases} \]

With Threshold Overlay

Various Thresholds

Threshold Histograms

Threshold Images

Segmenting Cells

Cell Colony plot of chunk unnamed-chunk-16

  • We can peform the same sort of analysis with this image of cells
  • This time we can derive the model from the basic physics of the system
    • The field is illuminated by white light of nearly uniform brightness
    • Cells absorb light causing darker regions to appear in the image
    • Lighter regions have no cells
    • Darker regions have cells

Different Threshold Values

Cell Colony

plot of chunk unnamed-chunk-18

Other Image Types

While scalar images are easiest, it is possible for any type of image \[ I(x,y) = \vec{f}(x,y) \]

The intensity can be described with two seperate or a single joint probability density function \[ P_{\vec{f}\cdot \vec{i}}(x,y), P_{\vec{f}\cdot \vec{j}}(x,y) \] With Threshold Overlay

Applying a threshold

A threshold is now more difficult to apply since there are now two distinct variables to deal with. The standard approach can be applied to both \[ I(x,y) = \begin{cases} 1, & \vec{f}_x(x,y) \geq0.25 \text{ and}\\ & \vec{f}_y(x,y) \geq0.25 \\ 0, & \text{otherwise} \end{cases} \]

This can also be shown on the joint probability distribution as With Threshold Overlay

Applying a threshold

Given the presence of two variables; however, more advanced approaches can also be investigated. For example we can keep only components parallel to the x axis by using the dot product. \[ I(x,y) = \begin{cases} 1, & |\vec{f}(x,y)\cdot \vec{i}| = 1 \\ 0, & \text{otherwise} \end{cases} \]

Looking at Orientations

We can tune the angular acceptance by using the fact \[ \vec{x}\cdot\vec{y}=|\vec{x}| |\vec{y}| \cos(\theta_{x\rightarrow y}) \] \[ I(x,y) = \begin{cases} 1, & \cos^{-1}(\vec{f}(x,y)\cdot \vec{i}) \leq \theta^{\circ} \\ 0, & \text{otherwise} \end{cases} \]

A Machine Learning Approach to Image Processing

Segmentation and all the steps leading up to it are really a specialized type of learning problem.

Returning to the ring image we had before, we start with our knowledge or ground-truth of the ring

plot of chunk unnamed-chunk-25

Which we want to identify the from the following image by using image processing tools

plot of chunk unnamed-chunk-26

What does identify mean?

  • Classify the pixels in the ring as Foreground
  • Classify the pixels outside of the ring as Background

How do we quantify this?

  • True Positive values in the ring that are classified as Foreground
  • True Negative values outside the ring that are classified as Background
  • False Positive values outside the ring that are classified as Foreground
  • False Negative values in the ring that are classified as Background

plot of chunk unnamed-chunk-27

We can then apply a threshold to the image to determine the number of points in each category

plot of chunk unnamed-chunk-28

Ring Threshold Example

Try a number of different threshold values on the image and compare them to the original classification

plot of chunk unnamed-chunk-30

Thresh TP TN FP FN
0.0 224 0 217 0
0.2 224 20 197 0
0.4 217 99 118 7
0.6 146 172 45 78
0.8 54 210 7 170
1.0 0 217 0 224

plot of chunk unnamed-chunk-32

Apply Precision and Recall

  • Recall (sensitivity)= \( TP/(TP+FN) \)
  • Precision = \( TP/(TP+FP) \)
Thresh TP TN FP FN Recall Precision
0.30 224 53 164 0 100 58
0.38 219 92 125 5 98 64
0.46 200 122 95 24 89 68
0.54 177 157 60 47 79 75
0.62 138 181 36 86 62 79
0.70 112 198 19 112 50 85

ROC Curve

Reciever Operating Characteristic (first developed for WW2 soldiers detecting objects in battlefields using radar). The ideal is the top-right (identify everything and miss nothing)

plot of chunk unnamed-chunk-34

Comparing Different Filters

We can then use this ROC curve to compare different filters (or even entire workflows), if the area is higher the approach is better.

plot of chunk unnamed-chunk-35

Different approaches can be compared by area under the curve

plot of chunk unnamed-chunk-36

True Positive Rate and False Positive Rate

Another way of showing the ROC curve (more common for machine learning rather than medical diagnosis) is using the True positive rate and False positive rate

  • True Positive Rate (recall)= \( TP/(TP+FN) \)
  • False Positive Rate = \( FP/(FP+TN) \)

These show very similar information with the major difference being the goal is to be in the upper left-hand corner. Additionally random guesses can be shown as the slope 1 line. Therefore for a system to be useful it must lie above the random line.

plot of chunk unnamed-chunk-37

Evaluating Models

Practical Example: Self-Driving Cars, Finding the Road

A more complicated example is finding the road in an image with cars and buildings

Tuk Tuk

From these images, an expert labeled the calcifications by hand, so we have ground truth data on where they are:

plot of chunk unnamed-chunk-39

Applying a threshold

We can perform the same analysis on an image like this one, again using a simple threshold to evalulate how accurately we identify the stripes. Since they are white, we can selectively keep the white components.

plot of chunk unnamed-chunk-40

plot of chunk unnamed-chunk-41

Examining the ROC Curve

Thresh TP TN FP FN Recall Precision
5 46933 37707 173131 149 100 21
10 46561 73225 137613 521 99 25
18 41374 105069 105769 5708 88 28
26 30625 131072 79766 16457 65 28
35 15866 153080 57758 31216 34 22
49 2580 176604 34234 44502 5 7

plot of chunk unnamed-chunk-43

plot of chunk unnamed-chunk-44

ROC: True Positive Rate vs False Positive

plot of chunk unnamed-chunk-45

Using other pieces of information

Since we know the stripes are on the road, we can add another criteria (just the lower half), we can improve the ROC curve substantially.

plot of chunk unnamed-chunk-47

plot of chunk unnamed-chunk-49

The resulting ROC curves

plot of chunk unnamed-chunk-50

Other Image Types

Going beyond vector the possibilities for images are limitless. The following shows a tensor plot with the tensor represented as an ellipse. \[ I(x,y) = \hat{f}(x,y) \]

Once the variable count is above 2, individual density functions and a series of cross plots are easier to interpret than some multidimensional density hypervolume.

Variable distributions

With Threshold Overlay

Multiple Phases: Segmenting Shale

  • Here we have a shale sample measured with X-ray tomography with three different phases inside (clay, rock, and air).
  • The model is that because the chemical composition and density of each phase is different they will absorb different amounts of x-rays and appear as different brightnesses in the image Shale Sample

Ideally we would derive 3 values for the thresholds based on a model for the composition of each phase and how much it absorbs, but that is not always possible or practical.

  • While there are 3 phases clearly visible in the image, the histogram is less telling (even after being re-scaled). plot of chunk unnamed-chunk-55

Multiple Segmentations

For this exercise we choose arbitrarily 3 ranges for the different phases and perform visual inspection plot of chunk unnamed-chunk-56

The relation can explicitly be written out as \[ I(x) = \begin{cases} \text{Void}, & 0 \leq x \leq 0.2 \\ \text{Clay}, & 0.2 < x \leq 0.5 \\ \text{Rock}, & 0.5 < x \end{cases} \] Segmented Images

Implementation

The implementations of basic thresholds and segmentations is very easy since it is a unary operation of a single image \[ f(I(\vec{x})) \] In mathematical terms this is called a map and since it does not require information from neighboring voxels or images it can be calculated for each point independently (parallel). Filters on the other hand almost always depend on neighboring voxels and thus the calculations are not as easy to seperate.

Implementation Code

Matlab

The simplist is a single threshold in Matlab:

thresh_img = gray_img > thresh

A more complicated threshold:

thresh_img = (gray_img > thresh_a) & (gray_img > thresh_b)

Python (numpy)

thresh_img = gray_img > thresh

Python

thresh_img = map(lambda gray_val: gray_val>thresh,
                gray_img)

Java

boolean[] thresh_img = new boolean[x_size*y_size*z_size];
for(int x=x_min;x<x_max;x++)
  for(int y=y_min;y<y_max;y++)
    for(int z=z_min;z<z_max;z++) {
      int offset=(z*y_size+y)*x_size+x;
      thresh_img[offset]=gray_img[offset]>thresh;
    }

In C/C++

bool* thresh_img = malloc(x_size*y_size*z_size * sizeof (bool));

for(int x=x_min;x<x_max;x++)
  for(int y=y_min;y<y_max;y++)
    for(int z=z_min;z<z_max;z++) {
      int offset=(z*y_size+y)*x_size+x;
      thresh_img[offset]=gray_img[offset]>thresh;
    }

Pitfalls with Segmentation

Partial Volume Effect

  • The partial volume effect is the name for the effect of discretization on the image into pixels or voxels.
  • Surfaces are complicated, voxels are simple boxes which make poor representations
  • Many voxels are only partially filled, but only the voxels on the surface
  • Removing the first layer alleviates issue

Partial Volume Effect

  • x -> val Partial Volume Effect

  • Shown as a function of radius

Radius Mean Intensity Sd Intensity
2.0 0.0311250 0.1623387
3.5 0.0960250 0.2770132
5.0 0.1956250 0.3825768
6.5 0.3311250 0.4510394
8.0 0.5017250 0.4819382
9.5 0.7072841 0.4258651

Morphology

Returning to the original image of a cross Cross With Threshold Overlay We can now utilize information from neighborhood voxels to improve the results. These steps are called morphological operations.

Like filtering the assumption behind morphological operations are

  • nearby voxels in real images are related / strongly correlated with one another
  • noise and imaging artifacts are less spatially correlated.

Therefore these imaging problems can be alleviated by adjusting the balance between local and neighborhood values.

Fundamentals: Neighborhood

A neighborhood consists of the pixels or voxels which are of sufficient proximity to a given point. There are a number of possible definitions which largely affect the result when it is invoked.

  • A large neighborhood performs operations over larger areas / volumes
    • Computationally intensive
    • Can smooth out features
  • A small neighborhood performs operations over small areas / volumes
    • Computationally cheaper
    • Struggles with large noise / filling large holes

The neighborhood is important for a large number of image and other (communication, mapping, networking) processing operations:

  • filtering
  • morphological operations
  • component labeling
  • distance maps
  • image correlation based tracking methods

Fundamentals: Neighbors in 2D

For standard image operations there are two definitions of neighborhood. The 4 and 8 adjacent neighbors shown below. Given the blue pixel in the center the red are the 4-adjacent and the red and green make up the 8 adjacent.

Formal Neighborhood Definitions

Formally these have two effectively equivalent, but different definitions.

  • Pixels which share a face (red line) with the pixel
  • Pixels which share a vertex (yellow dot) with the pixel

  • Pixels which are distance 1 from the pixel
  • Pixels which are distance \( \sqrt{2} \) from the pixel

Formal Neighborhood Definitions in 3D

In 3D there is now an additional group to consider because of the extra-dimension

  • Voxels which share a face (red line) with the voxel (6-adjacent)
  • Voxels which share an edge (yellow dot) with the voxel (18-adjacent)
  • Voxels which share a vertex (purple dot) with the voxel (26-adjacent)

  • Voxels which are distance 1 from the voxel
  • Voxels which are distance \( \sqrt{2} \) from the voxel
  • Voxels which are distance \( \sqrt{3} \) from the voxel

Erosion and Dilation

Erosion

If any of the voxels in the neighborhood are 0/false than the voxel will be set to 0

  • Has the effect of peeling the surface layer off of an object

Dilation

If any of the voxels in the neigbhorhood are 1/true then the voxel will be set to 1

  • Has the effect of adding a layer onto an object (dunking an strawberry in chocolate, adding a coat of paint to a car)

Applied Erosion and Dilation

Dilation

Taking an image of the cross at a too-high threshold, we show how dilation can be used to recover some of the missing pixels Dilation Example

Error in eval(expr, envir, enclos) : 
  could not find function "imgStdBinaryDilation"